install.packages("ggplot2")
install.packages("hexbin")
install.packages("foreign")
install.packages("grDevices")
install.packages("gmodels")
install.packages("lattice")
install.packages("arm")
install.packages("apsrtable")
install.packages("car")
install.packages("Zelig")
install.packages("lmtest")
install.packages("psych")
install.packages("effects")
install.packages("survival")
install.packages("RColorBrewer")
install.packages("plm")
install.packages("Ecdat")
install.packages("memisc")


############# Brian Forester ###################
############# MLE ##############################
############# Rage Replication #################

# Libraries and working directory
library(hexbin)
library(ggplot2) #load ggplot
library(foreign)
library(grDevices)
library(gmodels)
library(xtable)
library(lattice)
library(arm)
library(apsrtable)
library(car)
library(Zelig)
library(sandwich)
library(lmtest)
library(psych)
library(effects)
library(survival)
library(separationplot)
library(ROCR)
library(RColorBrewer)
library(plm)
library(Ecdat)
library(rms)
library(memisc)
library(aod)


setwd("C:/Users/bgforester/Desktop/Brian/2nd Semester/MLE/Rage_Against_Machines")  #Home
setwd("F:/2nd Semester/MLE/Rage_Against_Machines")  #Bunch
rm(list=ls())

data<-read.dta("Rage_RepData.dta")
head(data)
attach(data)

#Create table of summary statistics

sum.stats<-rbind(summary(wdl), summary(modern), summary(railway), summary(mech), 
                 summary(heli), summary(regime), summary(ltradegdp), summary(support), 
                 summary(power), summary(lenerpc), summary(ipatron), summary(lelev), 
                 summary(ldis),summary(numlang), summary(coldwar))

detach(data)

sd(na.omit(numlang))

colnames(sum.stats)<-c("Observations", "Standard deviation", "Median", "Mean", "Minimum", "Maximum", "Definition")
rownames(sum.stats)<-c("Outcome", "Modern", "Railway", "Mech", "Heli", "Regime", "Trade", "Support", "Power", "Energy", 
                       "Occupy", "Elevation", "Distance", "Language", "Cold war")

xtable(sum.stats)

#Investigate the dependent variable

modern.table<-table(data$modern,data$wdl)    
barplot(modern.table,  ylim=c(0,160),xlab="Outcome", ylab="Count", names.arg=c("Lose", "Tie", "Win"), 
        col=grey(c(0.25, 0.75)), las=1)
title("Counterinsurgency Outcomes")
legend("topleft",fill=grey(c(0.25, 0.75)), c("Not Modern (1800-1917)", "Modern (1918-2005)"), cex=0.75)


##########################################################################################
################################ REPLICATION OF RESULTS ##################################
##########################################################################################


#### Replication of Figure 1 ##############

data$quarter<-data$yearbg

data$quarter[data$quarter>=1800&data$quarter<=1825]<-LETTERS[1]
data$quarter[data$quarter>=1826&data$quarter<=1850]<-LETTERS[2]
data$quarter[data$quarter>=1851&data$quarter<=1875]<-LETTERS[3]
data$quarter[data$quarter>=1876&data$quarter<=1900]<-LETTERS[4]
data$quarter[data$quarter>=1901&data$quarter<=1925]<-LETTERS[5]
data$quarter[data$quarter>=1926&data$quarter<=1950]<-LETTERS[6]
data$quarter[data$quarter>=1951&data$quarter<=1975]<-LETTERS[7]
data$quarter[data$quarter>=1976&data$quarter<=2005]<-LETTERS[8]


win.table<-table(data$win, data$quarter)
win.percent<-prop.table(win.table,2)

names<-c("1800-\n1825", "1826-\n1850", "1851-\n1875", "1876-\n1900", "1901-\n1925",
         "1926-\n1950", "1951-\n1975", "1976-\n2005")

barplot(100*win.percent[2,], xlim=c(0.5,10), width=.8, space=.6, ylim=c(0,100), las=1, axes=F,
        names.arg=names, col="black", cex.names=0.8)

axis(side = 1, cex.axis=0.8, tck=0, labels=F)
axis(side = 2, at = seq(0, 100, 10), las=1,tck = -0.015, cex.axis=0.8)
box()

title(xlab="Years", ylab="Percentage won", cex.lab=0.8, line=2)


#### Replication of Figure 2 ##############


data$sold.per.vec<-c((1000*data$milper)/data$vec)

par(mar=c(3, 4, 4, 2)+0.1)

layout(rbind(1, 2), heights=c(7, 1)) 

plot(2005, 3000, type="n", xlim=c(1920,2005), ylim=c(0,3000), xlab="", ylab="",axes=F)

axis(side = 1, at = seq(1920, 2000, 20), tck = -0.015,cex.axis=0.8)
axis(side = 2, at = seq(0, 3000, 500), las=1,tck = -0.015, cex.axis=0.8)
axis(side = 3, at = seq(1920, 2000, 20), tck = 0, labels = FALSE)
axis(side = 4, at = seq(0, 3000, 500), tck = 0, labels = FALSE)
box()

points(data$yearbg, data$sold.per.vec, pch=20, cex=1.5)
lines(c(1920,2002), c(925, 110), lwd=1)

title(ylab="Soldiers per motor vehicle", cex.lab=0.8)
title(xlab="Year", line=2, cex.lab=0.8)

par(mar=c(0, 4, 0, 2)) 
plot.new() 

legend("topleft", pch=20, pt.cex=1, c("Incumbent mechanization level"), bty="n",cex=0.8)
legend("topright", lty=1,  c("Fitted values"), bty="n",cex=0.8)

layout(matrix(1))

par(mar=c(5, 4, 4, 2)+0.1)


########### Replication of the Models ####################
#Two fits per model - one without clustered SEs, and one with clustered SEs.
#I use the lrm function, from the "rms" package to compute the robust/clustered standard errors for 
#each model.  The clustered models are used for the remainder of the analysis.

#Model 1

model1<-polr(as.ordered(wdl)~modern, data=data, method="logistic", Hess=T)
summary(model1)

model1.r<-lrm(as.ordered(wdl)~modern, data=data, method="lrm.fit", x=T,y=T)
#Convert to robust standard errors
model1.r<-robcov(model1.r, data$ccode,method="huber")
print(model1.r)



#Model 2

model2<-polr(as.ordered(wdl)~modern+regime+support+power+lenerpc+occ+lelev+ldis+coldwar, 
            data=data, method="logistic", Hess=T)
summary(model2)

#Robust/Clustered SEs
model2.r<-lrm(wdl~modern+regime+support+power+lenerpc+occ+lelev+ldis+coldwar, 
             data=data, method="lrm.fit", x=T, y=T)
model2.r<-robcov(model2.r, data$ccode,method="huber")
print(model2.r)



#Model 3 (add railway for foraging era)

model3<-polr(as.ordered(wdl)~railway+regime+support+power+lenerpc+occ+lelev+ldis, 
             data=data[data$modern==0,], method="logistic", Hess=TRUE)
summary(model3)

#Robust/Clustered SEs
model3.r<-lrm(as.ordered(wdl)~railway+regime+support+power+lenerpc+occ+lelev+ldis, 
             data=data[data$modern==0,], method="lrm.fit", x=T, y=T)
model3.r<-robcov(model3.r, data$ccode[data$modern==0],method="huber")
print(model3.r)

#Model 4 (add railway for foraging era)

model4<-polr(as.ordered(wdl)~regime+support+power+lenerpc+occ+lelev+ldis+coldwar, 
             data=data, method="logistic", Hess=TRUE)
summary(model4)

#Robust/clustered SEs
model4.r<-lrm(as.ordered(wdl)~regime+support+power+lenerpc+occ+lelev+ldis+coldwar, 
             data=data[data$modern==1,], method="lrm.fit", x=T,y=T)
model4.r<-robcov(model4.r, data$ccode[data$modern==1],method="huber")
print(model4.r)

#Model 5 (mech only)

model5<-polr(as.ordered(wdl)~mech, 
             data=data[data$modern==1,], method="logistic", Hess=TRUE)
summary(model5)

#Robust/clustered SEs
model5.r<-lrm(as.ordered(wdl)~mech, 
             data=data[data$modern==1,], method="lrm.fit", x=T, y=T)
model5.r<-robcov(model5.r, data$ccode,method="huber")
print(model5.r)


#Model 6 (mech only)

model6<-polr(as.ordered(wdl)~mech+regime+support+power+lenerpc+occ+lelev+ldis, 
             data=data[120:286,], method="logistic", Hess=TRUE)
summary(model6)

#Robust/clustered SEs
model6.r<-lrm(as.ordered(wdl)~mech+regime+support+power+lenerpc+occ+lelev+ldis, 
             data=data[120:286,], method="lrm.fit", x=T, y=T)
model6.r<-robcov(model6.r, data$ccode[120:286], method="huber")
print(model6.r)


#Model 7 (post WW2)

model7<-polr(as.ordered(wdl)~mech+regime+support+power+lenerpc+occ+lelev+ldis+numlang+ltradegdp, 
             data=data[152:286,], method="logistic", Hess=TRUE)
summary(model7)

#Robust/clustered SEs
model7.r<-lrm(as.ordered(wdl)~mech+regime+support+power+lenerpc+occ+lelev+ldis+numlang+ltradegdp, 
             data=data[152:286,], method="lrm.fit", x=T, y=T)
model7.r<-robcov(model7.r, data$ccode[152:286],method="huber")
print(model7.r)


#Model 8 (heli only)

model8<-polr(as.ordered(wdl)~heli+regime+support+power+lenerpc+occ+lelev+ldis+numlang+ltradegdp, 
             data=data[152:286,], method="logistic", Hess=TRUE)
summary(model8)

#Robust/clustered SEs
model8.r<-lrm(formula=wdl~heli+regime+support+power+lenerpc+occ+lelev+ldis+numlang+ltradegdp, 
             data=data[152:286,], method="lrm.fit", x=T, y=T)
model8.r<-robcov(model8.r, data$ccode[152:286],method="huber")
print(model8.r)

wald.test(model3.r$var, model3.r$coefficients, Terms=1:10)

#### Create Table 2 #####

toLatex(mtable(model1, model2, model3, model4, digits=2), useBooktabs=F)

#### Create Table 3 #####

toLatex(mtable(model5, model6, model7, model8, digits=2), useBooktabs=F)

### Create coefficient plots for models 2, 6, and 8 to be used in body of article ###

#model 2
varnames<-c("", "Cutpoints","Modern", "Regime","Support", "Power", "Energy", "Occupy", 
            "Elevation", "Distance", "Cold War" )
coefplot(c(-model2.r$coefficients[1:2],model2.r$coefficients[3:11]), sqrt(diag(model2.r$var)),  
         varnames=varnames, main="", cex.var=0.6)
title(main="Model 2 \n (1800-2005)", cex.main=.8, line=1.5, cex.lab=0.6)
#pdf export size 4.5 x 5 in portrait

#model 6
varnames<-c("", "Cutpoints", "Mech", "Regime","Support", "Power", "Energy", "Occupy", 
            "Elevation", "Distance")
coefplot(c(-model6.r$coefficients[1:2],model6.r$coefficients[3:10]), sqrt(diag(model6.r$var)),  
         varnames=varnames, main="", cex.var=0.6)
title(main="Model 6 \n (1918-2005)", cex.main=.8, line=1.5)
#pdf export size 4.5 x 5 in portrait

#model 8
varnames<-c("","Cutpoints","Heli", "Regime","Support", "Power", "Energy", "Occupy", 
            "Elevation", "Distance", "Language", "Trade")
coefplot(c(-model8.r$coefficients[1:2],model8.r$coefficients[3:12]), sqrt(diag(model8.r$var)),  
         varnames=varnames, main="", cex.var=0.6)
title(main="Model 8 \n (1945-2005)", cex.main=.8, line=1.5)
#pdf export size 4.5 x 5 in portrait


########## Replication of Figure 3 (Model 6 conditional plot) ############

m6.data<-data[120:286,] #subset data to post-WW1 cases only


X<-cbind(m6.data$mech, m6.data$regime, m6.data$support, m6.data$power, 
         m6.data$lenerpc,m6.data$occ,m6.data$lelev, 
         m6.data$ldis)

pmat<-c()

for(i in 1:167){
  betas<-mvrnorm(n=10000, mu=c(model6$coefficients, model6$zeta), Sigma=vcov(model6))
  p<-matrix(NA, nrow=10000, ncol=2)
  
  for(j in 1:10000){
    
    p[j,1]<-1-plogis(betas[j,10]-X[i,]%*%betas[j,1:8])
    p[j,2]<-X[i,1]
    
  }
pmat<-rbind(pmat, p)
}


colnames(pmat)<-c("p.hat","mech")

layout(matrix(1)) #ensure plot parameters are correct

par(mar=c(5, 4, 4, 2)+0.1) #ensure plot parameters are correct

plot(4,.5, type="n", xlab="", ylab="", xlim=c(1,4), ylim=c(.2,.5), axes=F)
lines(spline(pmat[,2], pmat[,1], method="natural", ties=median))

axis(side = 1, at = seq(1, 4, 1), tck = 0.03,cex.axis=0.8)
axis(side = 2, at = seq(.2, .5, .05), las=1,tck = 0.03, cex.axis=0.8)
axis(side = 3, at = seq(1, 4, 1), tck = 0, labels = FALSE)
axis(side = 4, at = seq(.2, .5, .05), tck = 0, labels = FALSE)
box()

title(xlab="Mechanization level", cex.lab=0.8, line=2)
title(ylab="Probability of incumbent win", cex.lab=0.8)


##########################################################################################
################################ ASSESS MODEL FITS #######################################
##########################################################################################


########## Cross Validation of model 2 #############

set.seed(12)

sims<-1000

m2.cross<-data[sample(nrow(data), size=286),]
m2.cross$partitions<-c(rep(LETTERS[1:13],22))

expectedvaluemat <- c()

#mech at its median, occupy at median, all others at mean
scenario<-c(median(data$modern), mean(data$regime, na.rm=T), mean(data$support, na.rm=T), mean(data$power, na.rm=T), 
            mean(data$lenerpc, na.rm=T),median(data$occ),mean(data$lelev, na.rm=T), 
            mean(data$ldis, na.rm=T), median(data$coldwar))


for(i in 1:13){
  model<-polr(as.ordered(wdl)~modern+regime+support+power+lenerpc+occ+lelev+ldis+coldwar, 
              data=m2.cross[m2.cross$partitions!=LETTERS[i],], method="logistic")
  

head(betas)
  betas<-mvrnorm(n = sims, mu=c(model$coefficients, model$zeta),Sigma=vcov(model))
  
  expectedvalues<-matrix(rep(NA,1000*2), ncol=2)
  
  for(j in 1:sims){
    expectedvalues[j,1] <- 1-plogis(betas[j,11]-scenario%*%betas[j,1:9])
    expectedvalues[j,2] <- LETTERS[i]
  }
  
  expectedvaluemat <- rbind(expectedvaluemat, expectedvalues)
}

EVs <- data.frame(expectedvaluemat)
names(EVs) <- c("Exp.Value","Partition")
EVs <- transform(EVs, Exp.Value = as.numeric(as.character(Exp.Value)))


plot(1,2, xlab="", ylab="", axes=F, type="n", xlim=c(.2,.7), ylim=c(0, 6))

for(i in 1:13){
  lines(density(EVs$Exp.Value[EVs$Partition!=LETTERS[i]]), col=i)
}


axis(side = 1, at = seq(.2, .7, .1), tck = 0.03,cex.axis=0.8)
axis(side = 2, at = seq(0, 7, 1), las=1,tck = 0.03, cex.axis=0.8)
axis(side = 3, at = seq(.2, .7, .1), tck = 0, labels = FALSE)
axis(side = 4, at = seq(0, 7, 1), tck = 0, labels = FALSE)
box()

title(main="Model 2 Cross validation", cex.main=0.8,xlab="Probability of Win", 
      cex.lab=0.8, line=2)
title(ylab="Density", cex.lab=0.8, line=2)
legend("topright", col=c(1:13), title="Fold", lty=1, lwd=1, c(LETTERS[1:13]), cex=0.6)

########## Separation plot of model 2 ##############
## Models success of predicting incumbent wins ##

X<-cbind(data$modern, data$regime, data$support,data$power, 
         data$lenerpc, data$occ, data$lelev, 
         data$ldis, data$coldwar)
m2.p <- 1-plogis(model2$zeta[2]-X%*%model2$coefficients)

separationplot(as.vector(m2.p[2:286]), data$win[2:286], newplot=F)

###### ROC Plot of model 2 #######

pred <- prediction(m2.p[2:286], data$win[2:286])
perf <- performance(pred, "tpr","fpr")

plot(perf, colorize = FALSE, main=" Model 2 ROC Plot", cex.main=1, cex.xlab=.6, cex.ylab=.6, 
     yaxis.las=1)
abline(0,1)


########## Cross Validation of model 6 #############


set.seed(12)

sims<-1000

m6.cross<-m6.data[sample(nrow(m6.data), size=167),]
m6.cross$partitions<-c(rep(LETTERS[1:11],13), rep(LETTERS[12:13], 12))

expectedvaluemat <- c()

#mech at its median, occupy at median, all others at mean
scenario<-c(median(m6.data$mech), mean(m6.data$regime, na.rm=T), mean(m6.data$support, na.rm=T), mean(m6.data$power, na.rm=T), 
             mean(m6.data$lenerpc, na.rm=T),0,mean(m6.data$lelev, na.rm=T), 
             mean(m6.data$ldis, na.rm=T))


for(i in 1:13){
  model<-polr(as.ordered(wdl)~mech+regime+support+power+lenerpc+occ+lelev+ldis, 
              data=m6.cross[m6.cross$partitions!=LETTERS[i],], method="logistic", Hess=TRUE)

  betas<-mvrnorm(n = sims, mu=c(model$coefficients, model$zeta),Sigma=vcov(model))
  
  expectedvalues<-matrix(rep(NA,1000*2), ncol=2)
  
  for(j in 1:sims){
    expectedvalues[j,1] <- 1-plogis(betas[j,10]-scenario%*%betas[j,1:8])
    expectedvalues[j,2] <- LETTERS[i]
  }
  
  expectedvaluemat <- rbind(expectedvaluemat, expectedvalues)
}

EVs <- data.frame(expectedvaluemat)
names(EVs) <- c("Exp.Value","Partition")
EVs <- transform(EVs, Exp.Value = as.numeric(as.character(Exp.Value)))


plot(1,2, xlab="", ylab="", axes=F, type="n", xlim=c(.2,.75), ylim=c(0, 7))

for(i in 1:13){
lines(density(EVs$Exp.Value[EVs$Partition!=LETTERS[i]]), col=i)
}


axis(side = 1, at = seq(.2, .7, .1), tck = 0.03,cex.axis=0.8)
axis(side = 2, at = seq(0, 7, 1), las=1,tck = 0.03, cex.axis=0.8)
axis(side = 3, at = seq(.2, .7, .1), tck = 0, labels = FALSE)
axis(side = 4, at = seq(0, 7, 1), tck = 0, labels = FALSE)
box()

title(main="Model 6 Cross validation", cex.main=0.8,xlab="Probability of Win", 
      cex.lab=0.8, line=2)
title(ylab="Density", cex.lab=0.8, line=2)
legend("topright", col=c(1:13), title="Fold", lty=1, lwd=1, c(LETTERS[1:13]), cex=0.6)

########## Separation plot of model 6 ##############

X<-cbind(data$mech, data$regime, data$support,data$power, 
            data$lenerpc, data$occ, data$lelev, 
            data$ldis)
m6.p <- 1-plogis(model6$zeta[2]-X%*%model6$coefficients)

separationplot(as.vector(m6.p[120:286]), data$win[120:286], newplot=F)

###### ROC Plot of model 6 #######

pred <- prediction(m6.p[120:286], data$win[120:286])
perf <- performance(pred, "tpr","fpr")

plot(perf, colorize = FALSE, main=" Model 6 ROC Plot", cex.main=1, cex.xlab=.6, cex.ylab=.6, 
     yaxis.las=1)
abline(0,1)


########## Cross Validation of model 8 #############

m8.data<-data[152:286,]

m8.cross<-m8.data[sample(nrow(m8.data), size=135),]
m8.cross$partitions<-c(rep(LETTERS[1:9],15))

expectedvaluemat <- c()
m8.cross$partitions
#heli, occ at its median, occupy at median, all others at mean
scenario<-c(median(m8.data$heli, na.rm=T), mean(m8.data$regime, na.rm=T), mean(m8.data$support, na.rm=T), mean(m8.data$power, na.rm=T), 
            mean(m8.data$lenerpc, na.rm=T),median(m8.data$occ, na.rm=T),mean(m8.data$lelev, na.rm=T), 
            mean(m8.data$ldis, na.rm=T), mean(m8.data$numlang, na.rm=T), mean(m8.data$ltradegdp, na.rm=T))


for(i in 1:9){
  model<-polr(as.ordered(wdl)~heli+regime+support+power+lenerpc+occ+lelev+ldis+numlang+ltradegdp, 
              data=m6.cross[m6.cross$partitions!=LETTERS[i],], method="logistic", Hess=TRUE)
  
  betas<-mvrnorm(n = sims, mu=c(model$coefficients, model$zeta),Sigma=vcov(model))
  
  expectedvalues<-matrix(rep(NA,1000*2), ncol=2)
  
  for(j in 1:sims){
    expectedvalues[j,1] <- 1-plogis(betas[j,12]-scenario%*%betas[j,1:10])
    expectedvalues[j,2] <- LETTERS[i]
  }
  
  expectedvaluemat <- rbind(expectedvaluemat, expectedvalues)
}

EVs <- data.frame(expectedvaluemat)
names(EVs) <- c("Exp.Value","Partition")
EVs <- transform(EVs, Exp.Value = as.numeric(as.character(Exp.Value)))


plot(1,2, xlab="", ylab="", axes=F, type="n", xlim=c(.2,.75), ylim=c(0, 6))

for(i in 1:9){
  lines(density(EVs$Exp.Value[EVs$Partition!=LETTERS[i]]), col=i)
}


axis(side = 1, at = seq(.2, .7, .1), tck = 0.03,cex.axis=0.8)
axis(side = 2, at = seq(0, 6, 1), las=1,tck = 0.03, cex.axis=0.8)
axis(side = 3, at = seq(.2, .7, .1), tck = 0, labels = FALSE)
axis(side = 4, at = seq(0, 6, 1), tck = 0, labels = FALSE)
box()

title(main="Model 8 Cross validation", cex.main=0.8,xlab="Probability of Win", 
      cex.lab=0.8, line=2)
title(ylab="Density", cex.lab=0.8, line=2)
legend("topright", col=c(1:9), title="Fold", lty=1, lwd=1, c(LETTERS[1:9]), cex=0.6)


########## Separation plot of model 8 ##############

X<-cbind(data$heli, data$regime, data$support,data$power, 
         data$lenerpc, data$occ, data$lelev, 
         data$ldis, data$numlang, data$ltradegdp)
m8.p <- 1-plogis(model8$zeta[2]-X%*%model8$coefficients)

separationplot(as.vector(m8.p[152:286]), data$win[152:286], newplot=F)


###### ROC Plot of model 8 #######

pred <- prediction(m8.p[152:286], data$win[152:286])
perf <- performance(pred, "tpr","fpr")

plot(perf, colorize = FALSE, main=" Model 8 ROC Plot", cex.main=1, cex.xlab=.6, cex.ylab=.6, 
     yaxis.las=1)
abline(0,1)



###################################################################################
############################## EXTENSION (SIMS) ###################################
###################################################################################


########## Simulations using model 2 for modern/not modern ############


set.seed(12)

sims<-10000

pm<-rep(0,sims)
pnm<-rep(0,sims)

mean(pnh)

scenariom<-c(1, mean(data$regime, na.rm=T), mean(data$support, na.rm=T), mean(data$power, na.rm=T), 
             mean(data$lenerpc, na.rm=T),median(data$occ, na.rm=T),mean(data$lelev, na.rm=T), 
             mean(data$ldis, na.rm=T), median(data$coldwar))
scenarionm<-c(0, mean(data$regime, na.rm=T), mean(data$support, na.rm=T), mean(data$power, na.rm=T), 
             mean(data$lenerpc, na.rm=T),median(data$occ, na.rm=T),mean(data$lelev, na.rm=T), 
             mean(data$ldis, na.rm=T), median(data$coldwar))



betas<-mvrnorm(n = sims, mu=model2.r$coefficients, Sigma=vcov(model2.r))

for (i in 1:sims){
  pm[i]<-1-plogis(-betas[i,2]-scenariom%*%betas[i,3:11])
  pnm[i]<-1-plogis(-betas[i,2]-scenarionm%*%betas[i,3:11])

  
}

#Plot simulation results

plot(1,3,type="n", xlab="", ylab="", xlim=c(.1,1.0), ylim=c(-0.5,8), axes=F)

axis(side = 1, at = seq(0, 1, .1), tck = 0.03, cex.axis=0.8)
axis(side = 2, at = seq(0, 8, 1), las=1,tck = 0.03, cex.axis=0.8)
axis(side = 3, at = seq(0, 1, .1), tck = 0, labels = FALSE)
axis(side = 4, at = seq(0, 8, 1), tck = 0, labels = FALSE)
box()

lines(density(pm), col="#5E3C99", lwd=2) #modern
lines(density(pnm), col="#E66101", lwd=2)  #not modern

segments(quantile(pm, 0.025), -.2, quantile(pm, .975), -.2, lty=2, lwd=2, col="#998EC3")

segments(quantile(pnm, 0.025), -.35, quantile(pnm, .975), -.35, lty=2, lwd=2, col="#F1A340")

points(mean(pm), -.2, pch=20, col="#998EC3", cex=2)

points(mean(pnm), -.35, pch=20, col="#F1A340", cex=2)


title(main="Modern Simulation \n Model 2", cex.main=0.8)
title(xlab="Probability of Win", ylab="Density", cex.lab=0.8, line=2)
legend("topright",lwd=c(2,2,2), lty=c(1,1,2), col=c("#998EC3", "#F1A340","black"), 
       pch=c(NA,NA,20),c("machine era (1918-2005)", "foraging era (1800-1917)", "95% CI"), cex=0.6, pt.cex=1)



########## Simulations using model 6 for high/low levels of insurgent support ############


set.seed(12)

sims<-10000

p0<-rep(0,sims)
p1<-rep(0,sims)
p2<-rep(0,sims)



scenario0<-c(2, mean(m6.data$regime, na.rm=T), min(m6.data$support, na.rm=T), mean(m6.data$power, na.rm=T), 
             mean(m6.data$lenerpc, na.rm=T),median(m6.data$occ, na.rm=T),mean(m6.data$lelev, na.rm=T), 
             mean(m6.data$ldis, na.rm=T))
scenario1<-c(2, mean(m6.data$regime, na.rm=T), median(m6.data$support, na.rm=T), mean(m6.data$power, na.rm=T), 
             mean(m6.data$lenerpc, na.rm=T),median(m6.data$occ, na.rm=T),mean(m6.data$lelev, na.rm=T), 
             mean(m6.data$ldis, na.rm=T))
scenario2<-c(2, mean(m6.data$regime, na.rm=T), max(m6.data$support, na.rm=T), mean(m6.data$power, na.rm=T), 
             mean(m6.data$lenerpc, na.rm=T),median(m6.data$occ, na.rm=T),mean(m6.data$lelev, na.rm=T), 
             mean(m6.data$ldis, na.rm=T))



betas<-mvrnorm(n = sims, mu=model6.r$coefficients,Sigma=vcov(model6.r))

for (i in 1:sims){
  p0[i]<-1-plogis(-betas[i,2]-scenario0%*%betas[i,3:10])
  p1[i]<-1-plogis(-betas[i,2]-scenario1%*%betas[i,3:10])
  p2[i]<-1-plogis(-betas[i,2]-scenario2%*%betas[i,3:10])

}

#Plot simulation results

plot(1,3,type="n", xlab="", ylab="", xlim=c(0,0.85), ylim=c(-0.5,7.5), axes=F)

axis(side = 1, at = seq(0, 1, .1), tck = 0.03, cex.axis=0.8)
axis(side = 2, at = seq(0, 8, 1), las=1,tck = 0.03, cex.axis=0.8)
axis(side = 3, at = seq(0, 1, .1), tck = 0, labels = FALSE)
axis(side = 4, at = seq(0, 8, 1), tck = 0, labels = FALSE)
box()

lines(density(p0), col="#5E3C99", lwd=2) #Low support
lines(density(p2), col="#E66101", lwd=2)  #High support

segments(quantile(p0, 0.025), -.2, quantile(p0, .975), -.2, lty=2, lwd=2, col="#998EC3")

segments(quantile(p2, 0.025), -.35, quantile(p2, .975), -.35, lty=2, lwd=2, col="#F1A340")

points(mean(p0), -.2, pch=20, col="#998EC3", cex=2)

points(mean(p2), -.35, pch=20, col="#F1A340", cex=2)


title(main="Support Simulation \n Model 6", cex.main=0.8)
title(xlab="Probability of Win", ylab="Density", cex.lab=0.8, line=2)
legend("topright",lwd=c(2,2,2), lty=c(1,1,2), col=c("#998EC3", "#F1A340","black"), 
       pch=c(NA,NA,20),c("no insurgent support", "insurgent sanctuary and support", "95% CI"), cex=0.6, pt.cex=1)





########## Simulations using model 6 for high/low levels of mechanization ############


set.seed(12)

sims<-10000

p1<-rep(0,sims)
p2<-rep(0,sims)
p3<-rep(0,sims)
p4<-rep(0,sims)


scenario1<-c(1, mean(m6.data$regime, na.rm=T), mean(m6.data$support, na.rm=T), mean(m6.data$power, na.rm=T), 
            mean(m6.data$lenerpc, na.rm=T),0,mean(m6.data$lelev, na.rm=T), 
            mean(m6.data$ldis, na.rm=T))
scenario2<-c(2, mean(m6.data$regime, na.rm=T), mean(m6.data$support, na.rm=T), mean(m6.data$power, na.rm=T), 
            mean(m6.data$lenerpc, na.rm=T),0,mean(m6.data$lelev, na.rm=T), 
            mean(m6.data$ldis, na.rm=T))
scenario3<-c(3, mean(m6.data$regime, na.rm=T), mean(m6.data$support, na.rm=T), mean(m6.data$power, na.rm=T), 
            mean(m6.data$lenerpc, na.rm=T),0,mean(m6.data$lelev, na.rm=T), 
            mean(m6.data$ldis, na.rm=T))
scenario4<-c(4, mean(m6.data$regime, na.rm=T), mean(m6.data$support, na.rm=T), mean(m6.data$power, na.rm=T), 
            mean(m6.data$lenerpc, na.rm=T),0,mean(m6.data$lelev, na.rm=T), 
            mean(m6.data$ldis, na.rm=T))


betas<-mvrnorm(n = sims, mu=model6.r$coefficients,Sigma=vcov(model6.r))

for (i in 1:sims){
  p1[i]<-1-plogis(-betas[i,2]-scenario1%*%betas[i,3:10])
  p2[i]<-1-plogis(-betas[i,2]-scenario2%*%betas[i,3:10])
  p3[i]<-1-plogis(-betas[i,2]-scenario3%*%betas[i,3:10])
  p4[i]<-1-plogis(-betas[i,2]-scenario4%*%betas[i,3:10])
}

#Plot simulation results

plot(1,3,type="n", xlab="", ylab="", xlim=c(0.05,0.85), ylim=c(-0.5,7), axes=F)

axis(side = 1, at = seq(0, 1, .1), tck = 0.03, cex.axis=0.8)
axis(side = 2, at = seq(0, 7, 1), las=1,tck = 0.03, cex.axis=0.8)
axis(side = 3, at = seq(0, 1, .1), tck = 0, labels = FALSE)
axis(side = 4, at = seq(0, 7, 1), tck = 0, labels = FALSE)
box()

lines(density(p1), col="#5E3C99", lwd=2) #Low mech
lines(density(p4), col="#E66101", lwd=2)  #High mech

segments(quantile(p1, 0.025), -.2, quantile(p1, .975), -.2, lty=2, lwd=2, col="#998EC3")

segments(quantile(p4, 0.025), -.35, quantile(p4, .975), -.35, lty=2, lwd=2, col="#F1A340")

points(mean(p1), -.2, pch=20, col="#998EC3", cex=2)

points(mean(p4), -.35, pch=20, col="#F1A340", cex=2)


title(main="Mechanization Simulation \n Model 6", cex.main=0.8)
title(xlab="Probability of Win", ylab="Density", cex.lab=0.8, line=2)
legend("topright",lwd=c(2,2,2), lty=c(1,1,2), col=c("#998EC3", "#F1A340","black"), 
       pch=c(NA,NA,20),c("low mech", "high mech", "95% CI"), cex=0.6, pt.cex=1)



########## Simulations using model 8 for heli/no heli ############


m8.data<-data[152:286,]

sims<-10000

pnh<-rep(0,sims)
ph<-rep(0,sims)


scenarionh<-c(0, mean(m8.data$regime, na.rm=T), mean(m8.data$support, na.rm=T), mean(m8.data$power, na.rm=T), 
             mean(m8.data$lenerpc, na.rm=T),0,mean(m8.data$lelev, na.rm=T), 
             mean(m8.data$ldis, na.rm=T), mean(m8.data$numlang, na.rm=T), mean(m8.data$ltradegdp, na.rm=T))

scenarioh<-c(1, mean(m8.data$regime, na.rm=T), mean(m8.data$support, na.rm=T), mean(m8.data$power, na.rm=T), 
              mean(m8.data$lenerpc, na.rm=T),0,mean(m8.data$lelev, na.rm=T), 
              mean(m8.data$ldis, na.rm=T), mean(m8.data$numlang, na.rm=T), mean(m8.data$ltradegdp, na.rm=T))




betas<-mvrnorm(n = sims, mu=model8.r$coefficients,Sigma=vcov(model8.r))

for (i in 1:sims){
  pnh[i]<-1-plogis(-betas[i,2]-scenarionh%*%betas[i,3:12])
  ph[i]<-1-plogis(-betas[i,2]-scenarioh%*%betas[i,3:12])

}

#Plot simulation results

plot(1,3,type="n", xlab="", ylab="", xlim=c(0,.75), ylim=c(-0.5,8), axes=F)

axis(side = 1, at = seq(0, .75, .1), tck = 0.03, cex.axis=0.8)
axis(side = 2, at = seq(0, 8, 1), las=1,tck = 0.03, cex.axis=0.8)
axis(side = 3, at = seq(0, 8, .1), tck = 0, labels = FALSE)
axis(side = 4, at = seq(0, 6, 1), tck = 0, labels = FALSE)
box()

lines(density(pnh), col="#5E3C99", lwd=2) #No heli
lines(density(ph), col="#E66101", lwd=2)  #Heli

segments(quantile(pnh, 0.025), -.3, quantile(pnh, .975), -.3, lty=2, lwd=2, col="#998EC3")

segments(quantile(ph, 0.025), -.45, quantile(ph, .975), -.45, lty=2, lwd=2, col="#F1A340")

points(mean(pnh), -.3, pch=20, col="#998EC3", cex=2)

points(mean(ph), -.45, pch=20, col="#F1A340", cex=2)


title(main="Helicopter Simulation \n Model 8", cex.main=0.8)
title(xlab="Probability of Win", ylab="Density", cex.lab=0.8, line=2)
legend("topright",lwd=c(2,2,2), lty=c(1,1,2), col=c("#998EC3", "#F1A340","black"), 
       pch=c(NA,NA,20),c("no heli", "heli", "95% CI"), cex=0.6, pt.cex=1)



########## Simulations using model 6 for status as occupier ############

set.seed(12)

sims<-10000

pO<-rep(0,sims)
pNO<-rep(0,sims)


scenarioO<-c(2, mean(m6.data$regime, na.rm=T), mean(m6.data$support, na.rm=T), mean(m6.data$power, na.rm=T), 
             mean(m6.data$lenerpc, na.rm=T),1,mean(m6.data$lelev, na.rm=T), 
             mean(m6.data$ldis, na.rm=T))
scenarioNO<-c(2, mean(m6.data$regime, na.rm=T), mean(m6.data$support, na.rm=T), mean(m6.data$power, na.rm=T), 
              mean(m6.data$lenerpc, na.rm=T),0,mean(m6.data$lelev, na.rm=T), 
              mean(m6.data$ldis, na.rm=T))


betas<-mvrnorm(n = sims, mu=model6.r$coefficients,Sigma=vcov(model6.r))

for (i in 1:sims){
  pO[i]<-1-plogis(-betas[i,2]-scenarioO%*%betas[i,3:10])
  pNO[i]<-1-plogis(-betas[i,2]-scenarioNO%*%betas[i,3:10])
  
}

#Plot simulation results

plot(1,3,type="n", xlab="", ylab="", xlim=c(0,.70), ylim=c(-0.5,7.5), axes=F)

axis(side = 1, at = seq(0, 1, .1), tck = 0.03, cex.axis=0.8)
axis(side = 2, at = seq(0, 7, 1), las=1,tck = 0.03, cex.axis=0.8)
axis(side = 3, at = seq(0, 1, .1), tck = 0, labels = FALSE)
axis(side = 4, at = seq(0, 7, 1), tck = 0, labels = FALSE)
box()

lines(density(pO), col="#5E3C99", lwd=2) #Occupier
lines(density(pNO), col="#E66101", lwd=2)  #Non-Occupier

segments(quantile(pO, 0.025), -.2, quantile(pO, .975), -.2, lty=2, lwd=2, col="#998EC3")

segments(quantile(pNO, 0.025), -.35, quantile(pNO, .975), -.35, lty=2, lwd=2, col="#F1A340")

points(mean(pO), -.2, pch=20, col="#998EC3", cex=2)

points(mean(pNO), -.35, pch=20, col="#F1A340", cex=2)


title(main="Occupier Simulation \n Model 6", cex.main=0.8)
title(xlab="Probability of Win", ylab="Density", cex.lab=0.8, line=2)
legend("topright",lwd=c(2,2,2), lty=c(1,1,2), col=c("#998EC3", "#F1A340","black"), 
       pch=c(NA,NA,20),c("occupy", "non-occupy", "95% CI"), cex=0.6, pt.cex=1)



########## Simulations using model 6.1 for democracy/non-democracy ############
## I use the binary variable democ to differentiate between democracies and non-democracies

model61.r<-lrm(as.ordered(wdl)~mech+democ+support+power+lenerpc+occ+lelev+ldis, 
              data=data, method="lrm.fit", x=T, y=T)
model61.r<-robcov(model61.r, data$ccode, method="huber")
print(model61.r)

set.seed(12)

m61.data<-data[120:286,]

sims<-100000

pD<-rep(0,sims)
pND<-rep(0,sims)



scenarioD<-c(median(m61.data$mech),1, mean(m61.data$support, na.rm=T), mean(m61.data$power, na.rm=T), 
             mean(m61.data$lenerpc, na.rm=T),median(m61.data$occ),mean(m61.data$lelev, na.rm=T), 
             mean(m61.data$ldis, na.rm=T))
scenarioND<-c(median(m61.data$mech),0, mean(m61.data$support, na.rm=T), mean(m61.data$power, na.rm=T), 
              mean(m61.data$lenerpc, na.rm=T),median(m61.data$occ),mean(m61.data$lelev, na.rm=T), 
              mean(m61.data$ldis, na.rm=T))

betas<-mvrnorm(n = sims, mu=model61.r$coefficients,Sigma=vcov(model61.r))

for (i in 1:sims){
  pD[i]<-1-plogis(-betas[i,2]-scenarioD%*%betas[i,3:10]) #Democ
  pND[i]<-1-plogis(-betas[i,2]-scenarioND%*%betas[i,3:10]) #non-Democ
  
}



plot(1,3,type="n", xlab="", ylab="", xlim=c(.1,.75), ylim=c(-0.75,6.5), axes=F)

axis(side = 1, at = seq(.1, .75, .1), tck = 0.03, cex.axis=0.8)
axis(side = 2, at = seq(0, 6, 1), las=1,tck = 0.03, cex.axis=0.8)
axis(side = 3, at = seq(.1, .75, .1), tck = 0, labels = FALSE)
axis(side = 4, at = seq(0, 6, 1), tck = 0, labels = FALSE)
box()

lines(density(pD), col="#5E3C99", lwd=2) #democratic
lines(density(pND), col="#E66101", lwd=2)  #non-democratic

segments(quantile(pD, 0.025), -.3, quantile(pD, .975), -.3, lty=2, lwd=2, col="#998EC3")

segments(quantile(pND, 0.025), -.65, quantile(pND, .975), -.65, lty=2, lwd=2, col="#F1A340")

points(mean(pD), -.3, pch=20, col="#998EC3", cex=2)

points(mean(pND), -.65, pch=20, col="#F1A340", cex=2)


title(main="Democracy Simulation \n Model 6", cex.main=0.8)
title(xlab="Probability of Win", ylab="Density", cex.lab=0.8, line=2)
legend("topright",lwd=c(2,2,2), lty=c(1,1,2), col=c("#998EC3", "#F1A340","black"), 
       pch=c(NA,NA,20),c("democratic", "non-democratic", "95% CI"), cex=0.6, pt.cex=1)




#####################################################################################
######################## Extension using duration of conflict #######################
#####################################################################################

#Create the duration variable from differncing start date and end date

data$duration<-as.Date(data$enddate, "%d/%m/%Y")-as.Date(data$startdate, "%d/%m/%Y")
data$lduration<-log(as.numeric(data$duration))


m61.data<-data[120:286,]


####################### Duration interacted with Democracy using model 6.1 #####################
#re-specify model 6.1, including an interaction term between duration and democracy

model61.ra<-lrm(as.ordered(wdl)~mech+democ+support+power+lenerpc+occ+lelev+ldis+lduration+lduration*democ, 
               data=data, method="lrm.fit", x=T, y=T)
model61.ra<-robcov(model61.ra, data$ccode, method="huber")
print(model61.ra)

varnames<-c("", "Cutpoints","Mech", "Democ","Support", "Power", "Energy", "Occupy", 
            "Elevation", "Distance", "Duration", "Duration x \n Democ")

coefplot(c(-model61.ra$coefficients[1:2],model61.ra$coefficients[3:12]), sqrt(diag(model61.ra$var)),  
        varnames=varnames, main="", cex.var=0.55)

title(main="Model 6.1 \n (1918-2005)", cex.main=.8, line=1.5)


## Likelihood ratio test, to compare model 61.ra (duration) with model 6.r (without duration)

pchisq(-2*logLik(model6.r)+2*logLik(model61.ra), 2, lower.tail=F)


m61a.data<-data[c(-168,-193,-201),]
m61a.data<-m61a.data[120:283,] #(N=164, remove missing observations for duration)


########## Separation plot of model 6.1a ##############

X<-cbind(m61a.data$mech, m61a.data$democ, m61a.data$support,m61a.data$power, 
         m61a.data$lenerpc, m61a.data$occ, m61a.data$lelev, 
         m61a.data$ldis, m61a.data$lduration, m61a.data$lduration*m61a.data$democ)
m61a.p <- 1-plogis(-model61.ra$coefficients[2]-X%*%model61.ra$coefficients[3:12])

separationplot(as.vector(m61a.p), m61a.data$win, newplot=F)
#(use 3 x 7 in landscape)

###### ROC Plot of model 6.1a #######

pred <- prediction(m61a.p, m61a.data$win)
perf <- performance(pred, "tpr","fpr")

plot(perf, colorize = FALSE, main=" Model 6.1 ROC Plot", cex.main=1, cex.xlab=.6, cex.ylab=.6, 
     yaxis.las=1)
abline(0,1)


## Probability of winning as a function of duration, for both democ and non democ

X.d<-cbind(2, 1, mean(m61.data$support, na.rm=T), mean(m61.data$power, na.rm=T),
           mean(m61.data$lenerpc, na.rm=T), 0, mean(m61.data$lelev, na.rm=T), 
           mean(m61.data$ldis, na.rm=T), na.omit(m61.data$lduration), na.omit(m61.data$lduration))

set.seed(12)

sims<-10000

p.d<-matrix(NA, nrow=nrow(X.d), ncol=sims+4)

betas<-mvrnorm(n = sims, mu=model61.ra$coefficients,Sigma=vcov(model61.ra))

for(i in 1:nrow(X.d)){
  
  p.d[i,1]<-1-plogis(-model61.ra$coefficients[2]-X.d[i,]%*%model61.ra$coefficients[3:12])
  
  for(j in 1:nrow(betas)){
    
    p.d[i,j+1]<-1-plogis(-betas[j,2]-X.d[i,]%*%betas[j,3:12])
  }
  p.d[i,sims+2]<-quantile(p.d[i,2:sims], .05)
  p.d[i,sims+3]<-quantile(p.d[i,2:sims], .95)
}

p.d[,sims+4]<-X.d[,9]

p.d<-p.d[order(p.d[,sims+4]),]



### Non-Democracy ###

X.nd<-cbind(2, 0, mean(m61.data$support, na.rm=T), mean(m61.data$power, na.rm=T),
           mean(m61.data$lenerpc, na.rm=T), 0, mean(m61.data$lelev, na.rm=T), 
           mean(m61.data$ldis, na.rm=T), na.omit(m61.data$lduration), 0)

set.seed(12)

p.nd<-matrix(NA, nrow=nrow(X.nd), ncol=sims+4)


betas<-mvrnorm(n = sims, mu=model61.ra$coefficients,Sigma=vcov(model61.ra))

for(i in 1:nrow(X.nd)){
  
  p.nd[i,1]<-1-plogis(-model61.ra$coefficients[2]-X.nd[i,]%*%model61.ra$coefficients[3:12])
  
  for(j in 1:nrow(betas)){
    
    p.nd[i,j+1]<-1-plogis(-betas[j,2]-X.nd[i,]%*%betas[j,3:12])
  }
  p.nd[i,sims+2]<-quantile(p.nd[i,2:sims], .05)
  p.nd[i,sims+3]<-quantile(p.nd[i,2:sims], .95)
}

p.nd[,sims+4]<-X.nd[,9]

p.nd<-p.nd[order(p.nd[,sims+4]),]


## Plot log of duration and probability of winning


plot(9,1,type="n", xlab="", ylab="", xlim=c(2,10), ylim=c(0,1), axes=F)

axis(side = 1, at = seq(2, 10, 1), tck = 0.03, cex.axis=0.8)
axis(side = 2, at = seq(0, 1, .1), las=1,tck = 0.03, cex.axis=0.8)
axis(side = 3, at = seq(2, 10, 1), tck = 0, labels = FALSE)
axis(side = 4, at = seq(0, 1, .1), tck = 0, labels = FALSE)
box()

title(xlab="Logged duration of conflict (days)", cex.lab=0.8, line=2)
title(ylab="Pr(Win)",cex.lab=0.8, line=2.5)
title(main="Interaction of Democracy and Duration \n Model 6.1", cex.main=0.8, line=1)


polygon(x=c(p.d[,sims+4], rev(p.d[,sims+4])), y=c(p.d[,sims+2],rev(p.d[,sims+3])), 
        border=NA, col="#D8DAEB50")

lines(x=p.d[,sims+4], y=p.d[,1], lty=1, lwd=3, col="#5E3C99")

polygon(x=c(p.nd[,sims+4], rev(p.nd[,sims+4])), y=c(p.nd[,sims+2],rev(p.nd[,sims+3])), 
        border=NA, col="#F1A34050")

lines(x=p.nd[,sims+4], y=p.nd[,1], lty=1, lwd=3, col="#E66101")

legend("topright", lwd=3, col=c("#5E3C99","#E66101"), c("democracy", "non-democracy"), cex=0.6)

####save as pdf in landscape, with dimensions 6.5 x 7#####




